The silver blaze property for QCD with heavy quarks from the lattice 
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Lattice QCD at finite density suffers from a severe sign problem, which has so far prohibited 
simulations of the cold and dense regime. Here we study the onset of nuclear matter employing 
a three-dimensional effective theory derived by combined strong coupling and hopping expansions, 
which is valid for heavy but dynamical quarks and has a mild sign problem only. Its numerical 
evaluations agree between a standard Metropolis and complex Langevin algorithm, where the latter 
is free of the sign problem. Our continuum extrapolated data clearly show a first order phase 
transition building up at \ib ~uib as the temperature approaches zero. An excellent description of 
the data is achieved by an analytic solution in the strong coupling limit. 
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QCD at zero temperature is expected to exhibit the sil- 
ver blaze property: when a chemical potential for baryon 
number \xb is switched on, initially all observables should 
be completely independent of fis, despite the fact that 
the partition function explicitly depends on it. This 
changes abruptly once the chemical potential exceeds a 
critical value fig c , for which the baryon number jumps 
from zero to a finite value and a transition to a con- 
densed state of nuclear matter takes place. The reason 
for this behavior is the mass gap in the fermionic spec- 
trum, where the nucleon mass tub represents the low- 
est energy level that can be populated once /ib ~ tub- 
While this picture is easy to see in terms of energy lev- 
els of nucleons in a Hamiltonian language, it is elusive in 
the fundamental formulation of QCD thermodynamics in 
terms of a path integral. There, chemical potential en- 
ters through the Dirac operators of the quark fields, and 
hence all Dirac eigenvalues are shifted for any value of 
fiB- The silver blaze property thus requires some exact 
cancellations to take place as long as /is <mj. 

An analytic derivation of the silver blaze property from 
the path integral exists only for the related case of fi- 
nite isospin chemical potential fij = ji u = —/id [![, 
where Bose-Einstein condensation of pions sets in at 
fij = m,/2. A numerical demonstration of the effect 
by means of lattice QCD has so far been impossible due 
to the so-called sign problem. For finite baryon chemi- 
cal potential the action becomes complex, prohibiting its 
use in a Boltzmann factor for Monte Carlo approaches 
with importance sampling. Several approximate meth- 
ods have been devised to circumvent this problem. These 
are valid in the regime \i < T, where they give consistent 
results (for a recent review see 0). However, the cold 
and dense region of QCD has so far been inaccessible to 
lattice simulations. A method avoiding importance sam- 
pling is stochastic quantization, where expectation values 
are obtained from equilibrium distributions of stochastic 
processes governed by a Langevin equation 0]. While 



this works for several models with a sign problem [4|], 
it is not generally valid for complex actions Using 
Langevin dynamics, the silver blaze property has been 
numerically demonstrated for the Bose condensation of 
complex scalar fields Q- This was recently reproduced 
using a worm algorithm on the flux representation of the 
complex action, which is free of the sign problem [8| . 

In this work we show that cold and dense lattice QCD 
is accessible within a 3d effective theory of Polyakov 
loops, which has been derived from the full lattice theory 
with Wilson fermions by means of strong coupling and 
hopping parameter expansions 0, [Tc| ■ The pure gauge 
part reproduces the critical temperature T c of the decon- 
finement transition in the continuum limit to a few per- 
cent accuracy 0]. The theory was extended to include 
heavy but dynamical Wilson quarks. The sign problem 
of the resulting effective theory being under full control, 
the finite-temperature deconfincment transition, includ- 
ing its surface of endpoints, was located for all chemi- 
cal potentials. The critical quark mass corresponding to 
fiB = was again found to quantitatively agree with full 
4d Wilson simulations [lOj . The current restriction to 
large quark masses ensures the validity of the hopping 
expansion, we comment on possible extensions later. 

The lattice QCD partition function with Wilson gauge 
action S g [U] and / = 1, . . . , Nj quark flavours with Wil- 
son fermion matrix Q(nf,(j,f) can be written as 
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defining a 3d effective action by integration over the 
spatial link variables. The S'?' [W] depend on tempo- 
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ral Wilson lines W(x) = E[f=i ^for), and the Sf 
are Z(iV c )-symmetric while the S? are not. The cou- 
plings of the effective theory are functions of the tem- 
poral extent N T of the 4d lattice, the fundamental rep- 
resentation character coefficient u(j3) = /3/18 + 0((3 2 ) 
with lattice gauge coupling (3 = 2N c /g 2 and the hop- 
ping parameters k/, which for heavy quarks are related 
to the quark masses as Kf — exp(— amf)/2. Moreover, 
hif(/if) = hif(-fif). The couplings are then ordered by 
increasing powers of their leading contributions. Up to 
several non-trivial orders, the gauge sector is dominated 
by the nearest-neighbor interaction between Polyakov 
loops U = TrW(xi), 



-s; tf [W] 



JJ [1 + 2\KeUL*] 

<ij> 

X(u, N T >5) = u N - exp 



(2) 



N T I 4u 4 + I2u 5 - 14m 6 - 36m 7 



295 



-u° + 



1851 
W 



1055797 
5120 



(3) 



The convergence properties of the existing terms as well 
as explicit comparison with full 4d thermal simulations 
demonstrate that for j3 < 6.5 the pure gauge sector is un- 
der control and the effect of higher couplings is negligible 
for N T > 6 When fermions are present, /? is shifted 
by 0(k 4 ) corrections, which we neglect here. 

The Z(iV c )-breaking terms can be written as factors 
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Summing all windings of the temporal loops this reads 
Ai = \\dct[l + h lf Wtf [I + h lf w}] 2 ■ (5) 
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with the couplings 
hit = c f 
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C f = (2 Kf e^f) N ^ = e W-™/)/T 5 c f (fi f ) = C f (- H ). 
From now on we consider Nf — 1 and drop the index 
"/" (i.e. [i = /xg/3), which is sufficient to see the essential 
features. Finally we need meson and baryon masses, 



am M = -2 ln(2«) - 6k 2 - 24k 2 
arris = — 3 1ii(2k) — 18k z 
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Let us begin our analysis of the cold and dense regime 
with the combined static and strong coupling limit. In 



this case we have (3 = A = and the partition function 
factorizes into exactly solvable single site integrals: 



Z(f3 = 0) 
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The group integration only yields non-zero results if the 
trivial representation is contained in the products of 
loops. This results in the survival of hadronic degrees 
of freedom only, 
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We recognize the partition function of an ideal gas of 
baryons (~ C 3 ), mesons (CC) and composites of those, 
as already discussed in [llj. For finite chemical poten- 
tial and zero temperature, C —> and we are left with 
baryons, where we have reinstated N c to illustrate the 
meaning of the exponents. Prefactors are identified as 
spin degeneracy, i.e. we have a spin-3/2 quadruplet for 
the three-quark baryon and a spin zero baryon made of 
six quarks. 

The quark density is now easily calculated, 

T d 1 4N C C N ° + 2N C C 2N ^ , 

n= Vd-^ nZ= ^ l + AC^ + C' 2 ^ ■ (10) 

In the high density limit the expression reduces to 



lim (a 3 n) = 2N C = N c (a 3 n B ^ t ) 
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As required for fermions obeying the Pauli principle, the 
quark density in lattice units saturates once all available 
states per lattice site labeled by spin, color (and flavor 
for Nf > 1) are occupied. Note that summation over all 
windings of the Wilson lines is necessary in order to ob- 
tain a determinant in the form Eq. ([5]), while a truncation 
to finite order would not show saturation. Next, consider 
finite chemical potential and the zero temperature limit, 
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Thus the static strong coupling limit shows the silver 
blaze property, with zero quark density for /i < m and a 
jump to saturation density for fi > m, corresponding to 
a first order phase transition at quark chemical potential 
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fj, c = m. In the static strong coupling limit the baryon 
mass is arris = —3 ln(2«;) = 3am, i.e. the onset happens 
at fig — m B as required by the bounds in [12(. (Note 
that for static quarks, m-s/3 = m„/2). Moreover, in the 
dense phase the free energy scales with N c , consistent 
with the conjecture in jl3j . 

Since the step-function behavior is encoded in the ef- 
fective couplings hif, we expect it to be also present in 
the fully interacting theory at finite coupling and quark 
mass. However, in this case the partition function has 
to be computed numerically and for finite values of N T , 
i.e. the zero-temperature limit has to be approached nu- 
merically. The effective theory features a sign problem, 
which however is mild compared to that of the full 4d the- 
ory and can be overcome by reweighting methods using a 
standard Metropolis algorithm. For hif = the effective 
theory can be cast into a flux representation and simu- 
lated with the worm algorithm, without any sign prob- 
lem. The two approaches gave consistent results for the 
deconfinement transition at finite temperature and den- 
sity [ltj. Here we include /i 2 ^ 0, as it is the first coupling 
of quark-quark terms ~ L(x)L(y). In this case we could 
not find a flux representation free of the sign problem and 
instead employ complex Langevin simulations as an in- 
dependent check. Indeed, that algorithm has been shown 
to work for a simple SU (3)one-site model as well as QCD 
in the heavy dense limit [J], which have structures very 
similar to our effective theory. All our simulations sat- 
isfy the convergence criterion in terms of the Langevin 
operator specified in @. 

In order to reach continuum QCD, we work at small 
hopping parameters k<10 -3 , close to but not in the 
static limit. In this case we can use the non-perturbative 
beta-function of pure gauge theory for the lattice spacing 
in units of the Sommer parameter, <z(/3)/ro l4 |. Taking 
ro = 0.5 fm sets a physical scale for our lattices, while 
N T tunes temperature via T = (aN T )~ l . To begin, let 
us consider a temperature T — 10 McV and fix the pion 
mass to m w = 20 GeV. Because of the short Compton 
wave length, small lattices are sufficient. Once the lattice 
spacing is chosen, /3 is fixed and Eq. (J7J determines the 
corresponding «(/?). 

Fig. [1] shows the baryon density in lattice units as a 
function of chemical potential in units of the baryon mass 
for /3 = 5.7. It is consistent with zero until the chem- 
ical potential approaches tob/3, where a transition or 
crossover is clearly visible which quickly reaches satura- 
tion level. The rise in the baryon density is accompa- 
nied by a rise in the Polyakov loop. This feature was 
also seen in 4d Langevin simulations Q and in the chiral 
strong coupling limit in the staggered discretization . 
By contrast, in two-color QCD with lighter masses Bose 
condensation and the rise of the Polyakov loop appear to 
reflect two distinct transitions (ltjj . It is not clear to us 
whether the rise of the Polyakov loop signals deconfine- 
ment in the presence of matter. Evaluating (£*), (L) us- 
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FIG. 1: Baryon density ns/n_g, S at, Polyakov loop (L) 
and conjugate Polyakov loop (L*) as a function of [ib/tub 
obtained from Monte Carlo (N a = 3), complex Langevin 
(N s — 6) and the static strong coupling limit, respectively. 
Lattice parameters /3 = 5.7, « = 0.0000887, N T = 116 corre- 
spond to m M = 20 GeV, T = 10 MeV, a = 0.17 fm. 



ing Eq. ([5]) with (7 = 0, the Polyakov/conjugate loop gets 
screened by the third/second terms without changing the 
nature of the medium, which is hadronic. This also ex- 
plains why L* is screened before L when ii > 0, while the 
opposite happens for ii < 0. The ensuing decrease of the 
Polyakov and conjugate loops is a consequence of satura- 
tion: all color orientations get populated once the lattice 
approaches filling. All quantities in Fig. Q] agree between 
the Metropolis and the Langevin algorithms. The lat- 
ter becomes vastly superior once simulations on larger 
volumes are necessary. 

It is very striking that the numerical results are re- 
produced excellently by the analytic solution to the free, 
static hadron gas discussed earlier. That the static limit 
works well is easy to understand, since our quarks are 
exceedingly heavy and 0(k 2 ) corrections are tiny. What 
is less obvious is that a simulation at /3 = 5.7 is well 
approximated by the strong coupling limit, j3 = 0. The 
reason is that the effective coupling of the gauge sector 
A(j9 = 5.7, N T = 115) ~ 10~ 27 . This is an important ob- 
servation. The convergence of the strong coupling expan- 
sion depends on j3 and sufficient to allow for an accurate 
estimate of the convergence radius j3 c < 6 for N T < 16 
in However, for /3 in the same range, lowering tem- 
perature increases N T and thus improves convergence in 
two ways: we move away from the limiting convergence 
radius and u{0) < 1 gets suppressed by ever higher pow- 
ers. In other words, cold QCD is more amenable to the 
strong coupling expansion than hot QCD, and the pure 
gauge sector plays a negligible role for the dynamics. 

Simulations of the effective theory being cheap, we 
have computed the baryon density for nine gauge cou- 
plings 5.7 < (3 < 6.1, corresponding to lattice spacings 
0.17 fm > a > 0.07 fm. The scaling of the result in phys- 
ical units is shown in Fig. [5J Since the quark density is 
a derivative of the physical pressure with respect to an 
external parameter, it is a finite quantity that does not 
renormalize in a non-perturbative calculation. (The pres- 
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FIG. 2: Baryon density ns/m^ as a function of the lattice 
spacing a at T = 10 MeV and several values of the chemical 
potential. Crosses correspond to MC data (N s = 3), lines to 
the static strong coupling limit, Eq.©, (dashed) and includ- 
ing C(fi: 2 )-corrections (solid). 
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FIG. 3: Baryon density in the continuum. In the zero tem- 
perature limit a jump to nuclear matter builds up. 



sure requires subtraction of divergent vacuum energies, 
Pphys(T) = p(T)—p(0); however, these are /i- independent 
and d^pphys is finite.) Massive Wilson fermions have 0(a) 
lattice corrections, hence the continuum approach is 



ns,lat(A*) ns,cont(M) 



A{n)a + B(n)a 2 



(13) 



This behaviour is borne out by the data for a > 0.09 
fm in Fig. [2] On the other hand, for a — > we see 
a downward bend that violates scaling and signals that 
our truncated series in /3, k are no longer valid: as the 
lattice gets finer, j5 and k(/3) grow and our truncated 
effective action eventually must fail. Adding or removing 
0(k 2 ) corrections indeed affects this downward bend, but 
not the rest of the curve. Note that there is a trade-off 
between k and /3. The lighter we make the quark mass, 
the larger k for a given lattice spacing and the earlier the 
breakdown of the hopping expansion. Thus the scaling 
behavior of the baryon density tells us when our effective 
theory breaks down. 

For the very heavy quarks studied here, our series are 
stable and controlled for lattice spacings down to a > 0.09 
fm, which is just entering the regime with leading order 



lattice corrections. Cutting our data for a < 0.09-0.11 
fm, we perform continuum extrapolations based on five 
to seven lattice spacings by fitting to Eq. (fl"3|) . We have 
followed this procedure for four different temperatures, 
resulting in the number densities displayed in Fig. [3] 
Clearly, the silver blaze property is realized also in the 
interacting theory and away from the static limit as the 
zero temperature limit is approached. Interestingly, the 
density at onset, when expressed in units of m^, is of the 
same order of magnitude as the physical nuclear density 
- 0.16 far 3 w 0.15 ■ 10~ 2 m 3 ro ton- 

For light quarks our truncated hopping series is not 
reliable. Since the hopping expansion provides a conver- 
gent series, the discontinuity in the couplings hif should 
still imply a transition to all orders, its location being 
shifted by corrections. Whether the pure gauge sector is 
similarly suppressed in that case remains to be seen. We 
plan to address this issue together with K 4 -corrections, a 
thorough discussion of the Langevin simulations and the 
physical dynamics in a future publication. 
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